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Abstract 

A framework for constructing integral preserving numerical schemes 
for time-dependent partial differential equations on non-uniform grids 
is presented. The approach can be used with both finite difference 
and partition of unity methods, thereby including finite element meth¬ 
ods. The schemes are then extended to accommodate r-, h- and 
p-adaptivity. To illustrate the ideas, the method is applied to the 
Korteweg-de Vries equation and the sine-Gordon equation. Results 
from numerical experiments are presented. 


1 Introduction 

Courant, Friedrichs and Lewy introduced difference schemes with conserva¬ 
tion properties in [^, where a discrete conservation law for a finite difference 
approximation of the wave equation was derived. Their methods are often 
called energy methods or energy-conserving methods [^, although the 
conserved quantity is often not energy in the physical sense. The primary 
motivation for developing conservative methods was originally to devise a 
norm that could guarantee global stability. This was still an objective, in 
addition to proving existence and uniqueness of solutions, when the energy 
methods garnered newfound interest in the 1950s and 1960s, resulting in 
new developments such as generalizations of the methods and more differ¬ 
ence schemes, summarized by Richtmyer and Morton in [^. In the 1970s, 
the motivation behind studying schemes that preserve invariant quantities 
changed, as the focus shifted to the conservation property itself. Li and Vu- 
Quoc presented in a historical survey of conservative methods developed 
up to the early 1990s. They state that this line of work is motivated by the 
fact that in some situations, the success of a numerical solution will depend 
on its ability to preserve one or more of the invariant properties of the original 
differential equation. In addition, as noted in [^[^, there is the general idea 
that transferring more of the properties of the original continuous dynamical 
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system over to a discrete dynamical system may lead to a more accurate nu¬ 
merical approximation of the solution, especially over long time intervals. In 
recent years, there has been a greater interest in developing systematic tech¬ 
niques applicable to larger classes of differential equations. Hairer, Lubich 
and Wanner give in a presentation of geometric integrators for differential 
equations, i.e. methods for solving ordinary differential equations (ODEs) 
that preserve a geometric structure of the system. Examples of such geomet¬ 
ric structures are symplectic structures, symmetries, reversing symmetries, 
isospectrality. Lie group structure, orthonormality, first integrals, and other 
invariants, such as volume and invariant measure. 

In this paper we will be concerned with the preservation of first integrals 
of PDEs. From the ODE literature we find that the most general methods 
for preserving first integrals are tailored schemes, in the sense that the vec¬ 
tor field of the ODE does not by itself provide sufficient information, so the 
schemes make explicit use of the first integral. An obvious approach in this 
respect is projection, where the solution is first advanced using any consistent 
numerical scheme and then this approximation is projected onto the appro¬ 
priate level set of the invariant. In the same class of tailored methods one 
also has the discrete gradient methods, usually attributed to Gonzalez [^. 
For the subclass of canonical Hamiltonian systems, the energy can be pre¬ 
served by means of a general purpose method called the averaged vector field 
method, see e.g. [^. 

The notion of discrete gradient methods for ordinary differential equa¬ 
tions has a counterpart for partial differential equations called the discrete 
variational derivative method. Such schemes have been developed since the 
late 1990s in a number of articles by Japanese researchers such as Furihata, 
Matsuo, Sugihara, and Yaguchi. A relatively recent account of this work can 
be found in the monograph [^. More recently, the development of integral 
preserving schemes for PDEs has been systematised and eased, in particular 
by using the aforementioned tools from ordinary differential equations, see 
for instance |10||11| . Most of the schemes one finds in the literature are based 
on a finite difference approach, and usually on fixed, uniform grids. There are 
however some exceptions. Yaguchi, Matsuo and Sugihara presented in [I^13| 
two different discrete variational derivative methods on fixed, non-uniform 
grids, specifically defined for certain classes of PDEs. Non-uniform grids 
are of particular importance for multidimensional problems, since the use 
of uniform grids will greatly restrict the types of domains possible to dis¬ 
cretize. Another important consequence of being able to use non-uniform 
grids is that it allows for the use of time-adaptive spatial meshes for solv¬ 
ing partial differential equations. Adaptive energy preserving schemes for 
the Korteweg-de Vries and Cahn-Hilliard equations have been developed re¬ 
cently |14] by Miyatake and Matsuo. The main objective of this paper is to 
propose a general framework for numerical methods for PDEs that combine 
mesh adaptivity with first integral conservation. 
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Several forms of adaptive methods exist; they can roughly be categorized 
as r-, h- and p-adaptive. When applying r-adaptivity, one keeps the number 
of degrees of freedom constant while modifying the mesh at each time step 
to e.g. cluster in problematic areas such as boundary layers or to follow 
wave fronts. When applying the Finite Difference Method (FDM) or the 
Finite Element Method (FEM), moving mesh methods may be used for r- 
adaptivity, some examples of which may be found in [l5| - [l7| . When using 
Partition of Unity Methods (PUM) (and in particular when using FEM), h- 
and p-adaptivity relate to adjusting the number of elements and the basis 
functions used on the elements, respectively. For PUM methods there exist 
strategies for h- and p-adaptivity based both on a priori and a posteriori 
error analysis 18 . Common to all of these strategies is that, based on 


estimated function values in preceding time steps, one can suggest improved 
discretization parameters for the next time step. In the FDM approach, these 
discretization parameters consist of the mesh points x, while in the PUM 
approach the parameters encompass information about both the mesh and 
the basis functions. We will, in general, denote a collection of discretization 
parameters by p, and assume that the discretization parameters are changed 
separately from the degrees of freedom u of the problem when using adaptive 
methods. That is, starting with an initial set of discretization parameters 
p® and initial values u^, one would hrst decide upon p^ before calculating 
u^, then Ending p^, then u^, etc., in a decoupled fashion. 

A first integral of a PDE is a functional X on an infinite-dimensional 
function space, whereas the numerical methods considered here will reduce 
the problem to a hnite-dimensional setting. Therefore, we cannot preserve 
the exact value of the first integral; instead, we will preserve a consistent 


approximation to the hrst integral, Xp 


u 


The approximation will be de¬ 


pendent on the discretization parameters p and, since adaptivity alters the 
discretization parameters, we will therefore aim to preserve the value of the 
approximated hrst integral across all discretization parameters, i.e. we will 
require that Xpn+i = Xp^ (u”). Here, and in the following, superscripts 

denote time steps unless otherwise specihed. 

In this article, we will present a method for developing adaptive numerical 
schemes that conserve an approximated hrst integral. In Section 2, the PDE 
problem is stated, and two classes of hrst integral preserving methods using 
arbitrary yet constant discretization parameters are presented; one using an 
FDM approach and the other a PUM approach for spatial discretization. A 
connection to previously existing methods is then established. In Section 3, 
we present a way of adding adaptivity to the methods from Section 2 and the 
modifications needed to retain the hrst integral preservation property, before 
showing that certain projection methods form a subclass of the methods thus 
obtained. Section 4 contains examples of the application of the methods to 
two PDFs and numerical results pertaining to the quality of the numerical 
solutions as compared to a standard implicit method. 
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2 Spatial discretization with fixed mesh 


2.1 Problem statement 

Consider a partial differential equation 

ut = f{x,u-^), xeOC 


u £B L^, 


( 2 . 1 ) 


where denotes u itself and its partial derivatives of any order with respect 
to the spatial variables xi, ....,Xd- We shall not specify the space B further, 
but assume that it is sufficiently regular to allow all operations used in the 
following. For ease of reading, all t-dependence will be suppressed in the 
notation wherever it is irrelevant. Also, from here on, square brackets are 
used to denote dependence on a function and its partial derivatives of any 
order with respect to the independent variables t and xi, ...,Xd- We recall the 
definition of the variational derivative of a functional H[ii\ as the function 
^[u] satisfying 


8H \ d 
— [u\,v) = — 


6u 


L2 


de 


H[u + ev] Vu G B, 


( 2 . 2 ) 


e=0 


and define a first integral of (2.1) to be a functional X\u] satisfying 

= 0, Vtt G B. 


^[u],/(x,u‘^) 


6u 


L2 


We may observe that X\u] is preserved over time, since this implies 


dX n du\ 


dt 


Furthermore, we may observe that if there exists some operator ^(x, n'^), 
skew-symmetric with respect to the inner product, such that 

/(x,W^) = S{x,u-^)^[u], 


then X[n] is a first integral of (2.11, and we can state (2.1) in the form 

<5X, 


ut = ^(x, n-^) —[u]. 

ou 


(2.3) 


This can be considered as the PDF analogue of an ODE with a first integral, 
in which case we have a system 


(ix 

- = S(x)VI(x), 


(2.4) 


where S'(x) is a skew-symmetric matrix [^. Note that Hamiltonian equa¬ 
tions are contained of this class of ODEs. For such differential equations. 
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there exist numerical methods preserving the hrst integral I (x), for instance 
the discrete gradient methods, which are of the form 
^n+l _ „n _ _ 

- - -= 5(x-,x"+i)V/(x",x-+i), 

where ,S'(x"', is a consistent skew-symmetric time-discrete approxima¬ 

tion to S{x) and V/(v,u) is a discrete gradient of /(x), i.e. a function 
satisfying 


(V/(v,u))'^(u-v) = /(u)-/(v), 
V/(u,u) = V/(u). 


(2.5) 

( 2 . 6 ) 


There are several possible choices of discrete gradients available, one of which 
is the Average Vector Field (AVF) discrete gradient |10| , given by 

1 

V/(v,u) = I V/(eu + (l-e)v)d^, 

0 


which will be used for numerical experiments in the hnal chapter. Our 
approach to solving ( 2.1[ ) on non-uniform grids is based upon considering 
the PDF in the form (2.3), reducing it to a system of ODEs of the form (2.4) 
and applying a discrete gradient method. This is done by Ending a discrete 
approximation Ip to I and using this to obtain a discretization in the spatial 
variables, which is achieved through either a hnite difference approach or a 
variational approach. 


2.2 Finite difference method 


In the finite difference approach, we restrict ourselves to obtaining approx¬ 
imate values of u at the grid points xq, which can be interpreted as 

quadrature points with some associated nonzero quadrature weights kq; 

The grid points constitute the discretization parameters p. We can then ap¬ 
proximate the inner product by quadrature to arrive at a weighted inner 
product: 


{u,v)j^ 2 = / u(x)n(x)dx 


M 

E 

1=0 


KiU{Xi)v{Xi) = u^D{k)v = (u, v)^ , 


where D(k) = diag(/{o, km)- Assume that there exists a consistent ap¬ 
proximation 2p(u) to the functional I[u\, dependent on the values of u at the 
points Xj. Then, we can characterize the discretized variational derivative 
by asserting that 


iXp, , \ d 


Xp(u -I- ev) Vv G 


pM+l 


e=0 
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meaning 


\ 5u 


^(u) Z)(k)v= (VXp(u))'v VvG 


oM+l 


from which we conclude that 

51, 


^(u) = D(k) ^VXp(u). 


(2.7) 


Using this as a discretization of [tt] and approximating ^(x, w^) by a matrix 
S'd(u), skew-symmetric with respect to (•, •)«, we obtain a discretization of 


(2.3) as: 


^ = 5p(u)VXp(u), 


( 2 . 8 ) 


where Sp{u) = >S'^(u)D(ai:) This system of ODEs is of the form (2.4), 


since 


Sp{nf = {Sdiu)DiK)-Y 

= D{Kr^Sd{ufD{K)D{Kr^ 
= -D{K)-^D{K)Sd{vL)D{Kr^ 
= -Sd{n)D{K)-^ 

= -5'p(u). 


This allows us to apply first integral preserving methods for systems of ODEs 
to solve the spatially discretized system. For example, we may consider using 
a discrete gradient VXp, and a skew-symmetric, time-discrete approximation 
5p(u”', u"'^^) to 5'p(u), where u"" = u(t„), tn = nAt. Then, the following 
scheme will preserve the approximated first integral Xp in the sense that 
Xp(u-+i) =Xp(u-): 


U^+l _ 

At 


= 5p(u-,u-+i)VXp(u-,u 


n 


(2.9) 


2.3 Partition of unity method 


One may also approach the problem of spatially discretizing the PDE through 
the use of variational methods such as the Partition of Unity Method (PUM) 
|20| , which generalizes the Finite Element Method (FEM). Here, the varia¬ 
tional structure of the functional derivative can be utilized in a natural way, 
such that one avoids having to approximate ^(x, tf^). We begin by stating 


a weak form of (2.3). Then, the problem consists of finding u € B such that 


{ut,v)^2 = ((^5(x,?r‘^)^[n],: 


V ) = 

L2 




Vr; G B. 

( 2 . 10 ) 
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Employing a Galerkin formulation, we restrict the search to a finite dimen¬ 
sional subspace = span{(/9o, ■■■fM} ^ and approximate u by the func¬ 
tion 

M 

u^{x,t) = '^Ui{t)ipi{x). 
i=0 

We denote by p the collection of discretization parameters defining B^; this 
includes information about mesh points, element types and shapes of basis 
functions. Furthermore, we define the canonical mapping <hp : 1^-^+^ — B^ 
given by 


M 


^p(u) = 

i=0 


( 2 . 11 ) 


and the discrete first integral Ip by 

Xp(u)=X(cl>p(u)). 

The following lemma will prove useful later in the construction of the method: 
Lemma 1. For any u^,v G B^, 


de 


X(u^ + ev) = (VXp(u))^v. 


e=0 


Proof. 

de 


e=0 


X{yh + eu) = ^ 


X(<hp(u + ev)) 

e=0 

= (^|$p(u + ev)],i5.p(u + fv) 


de 


L2 


e=0 


= <^^[^p(u + ev)],(V^p(u + ev))^v^ 


L2 


e=0 


SI 

5u 


[<i>p(u)],(V<l>p(u))^v 


L2 


M 


^ .■_A ^ 


i=0 


i=0 


□ 
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We observe that for u,v G the L? inner product has a discrete counter¬ 
part: 


M M 

{U,V)^2 = = (u,v)^ 

i=0 j=0 

with the symmetric positive definite matrix A given by Aij = {cpi, cpj)j^ 2 - 
Note also that equation (2.10) is satisfied in B^ if it is satisfied for all basis 
functions (pj. The Galerkin form of the problem therefore consists of finding 
Ui{t) such that 

^ J j S'J' \ 

E ^ {“.P-12) 


i=0 


L2 


This weak form is rather unwieldy and does not give rise to a system of the 
form (2.4), so in order to make further progress, we consider the projection 
of |f[r?Tonto B^: 


61^ 


M 


M 




6u 

where rci(u) = r(;^[<l>(u)] = w^[u^\ are coefficie nts th at will be characterized 
later. Replacing by its projection in (2.12) gives the approximate 

weak form: 


2 = 0 2 = 0 

,,/2 L./ll 


M 


2=0 2=0 

Thus, we obtain a system of equations for the coefficients up. 


M 


A^ = -R(u)w(u), 


(2.13) 


with the skew-symmetric matrix i?(u) given by B{u.)ji = (^pi, ^(x, 
Furthermore, we may characterize the vector w(u) by the following argu¬ 
ment: 


'61^ 




5T, 


d 


I{u^ -F ev) = (VXp(u))^v, 


e=0 


where the last equality holds by Lemma This holds for all v G and 

thus 


w(u) = ^ ^VXp(u). 


(2.14) 
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Inserting (2.14) into (2.13) and left-multiplying by A we are left with an 
ODE for the coefficients n,-: 


^ = S’p(u)VXp(u). 


(2.15) 


Here, 5'p(u) = —is a skew-symmetric matrix, and the system 
is thereby of the form (2.4), meaning Zp can be preserved numerically using 
e.g. discrete gradient methods as in equation (2.9). 


2.4 Discrete variational derivative methods 


Let us now define a general framework for the discrete variational derivative 
methods that encompass the methods presented by Furihata, Matsuo and 
coauthors in a number of publications including [^[^|12[[I^|21]. 


Definition 1. Let Zp be a consistent approximation to the functional Z\v\ 
discretized on p given by grid points Xj and quadrature weights Hi, i = 
0, ...,M. Then 5 (v^(v,u) is a discrete variational derivative of Xp(u) if it 
is a continuous function satisfying 




(l(v,u)’ 


u — V 


<5Xn 


(5(u,u) 


Xp(u) -Xp(v), 

(2.16) 

/>X„ 


(2.17) 


and the discrete variational derivative methods for solving PDEs on the form 


(2.3) are given by 


u”+i - u" 
At 


= 5d(u-,u"+i 


<5Xn 




(2.18) 


where S'd(u"', u^^^) is a time-discrete approximation to ^^(u), and itself 
skew-symmetric with respect to the inner product (•,•)«. 


Proposition 1. A discrete gradient method (2.9) applied to the system of 
ODEs (2.8) or (2.15) is equivalent to a discrete variational derivative method 
as given by (2.18\), with 


5,(u-,u-+1) = 5p(u-,u"+1)D(k), 


and the discrete variational derivative 


6Zp 

(5(v,u) 


D(Ac)-iVXp(v,u) 


satisfying (2.16)-(2J^. 


(2.19) 
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Proof. Applying (2.5), we get that, for the discrete variational derivative 
dehned by (2.19), 


^VXp(v,u),u-v)^ 


= {D{k) ^VXp (v, u))'^ D{k) (u - v) 

= VXp (v, u)'^ (u - v) = Xp(u) - Xp(v) 


and hence (2.16) is satished. Furthermore, applying (2.6) and (2.7), 
= X)(k)“^VXp(u, u) = D{k)-^VIp (u) = ^ (u) 
and (2.17[) is also satished. 


□ 


Consequently, all discrete variational derivative methods as given by 


(2.18) can be expressed as discrete gradient methods on the system of ODEs 
(2.8) or (2.15) obtained by discretizing (2.3) in space, and vice versa. 


3 Adaptive discretization 

3.1 Mapping solutions between parameter sets 

Assuming that adaptive strategies are employed, one would obtain a new 
set of discretization parameters p at each time step. After such a p has 
been found, the solution using the previous parameters must be transferred 
to the new parameter set before advancing to the next time step. This 
transfer procedure can be done in either a preserving or a non-preserving 
manner. Let p”, u”, p”+^ and denote the discretization parameters 

and the numerical values obtained at the current time step and next time 
step, respectively. Also, let u denote the values of u” transferred onto p"'^^ 
by whatever means. We call the transfer operation preserving if Xpn+i(u) = 
Xpn(u’^). If the transfer is preserving, then the next time step can be taken 
with a preserving scheme, e.g. the scheme 

n+l _ _ 

= V+du,u"+i)VXp.+i(u,u-+i), 
which is preserving in the sense that 
Xp„ + l(u"+l)-Xpn(u^) =Xpn + l(u"+l) -Xp„ + l(u) 

= (VXpn+i(u,u"'+^),u"'+^-u) 

= At (VXprL + l(u, 5pn + l(u, u"'^^)VXprL + l(u, 

= 0 , 

since Spn+i (u, is skew-symmetric. If non-preserving transfer is used, 

corrections are needed in order to obtain a preserving numerical method. 
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Proposition 2. The scheme 


u'^+i = U- V(^ ))^, ^At5pn+1 (u, u"+i )VXpn+i (u, ) 


(VXpn + l(u, 0*^+1 ),z) 


(3.1) 


where z is an arbitrary vector, is first integral preserving in the sense that 
Xpn + l(u”+l)-Xpn(u^) =0. 

Proof. 

Xp„+i(u"+i)-Xpn(u")=Xp„+i(u"+i) -Xp„+i(u) +Xp„+i(u) -Xpn(u”) 

= (VXpn+i(u,u"'+^),u”+^-u) +Xpn+l(u) -Xpu(u”) 

\ P ^ (VXp.+i(u,u-+i),z)/ 

= 0 . 


□ 


The correcting direction z should be chosen so as to obtain a minimal correc¬ 
tion, and such that (VXpn+i(u, z) ^ 0. One possibility is simply tak¬ 
ing z = VXpn+i (u, In the FDM case one may alternatively choose z = 

X)(K)“^VXpn+i (u, and in the PUM case, z = ^“^VXpn+i (u, 

When using the PUM formulation, one may obtain a method for pre¬ 
serving transfer in the following manner. Any changes through e.g. r- p- 
and/or /i-refinement between time steps will result in a change in the shape 
and/or number of basis functions. Denote by = span{(/jj}(^Q the trial 
space from the current time step and by = span{(^j}(^Q the trial space 
for the next time step, and note that in general, M ^ M. We do not concern 
ourselves with how the new basis is found, but simply acknowledge that the 
basis changes through adaptivity measures as presented in e.g. 15 or 18 . 
Our task is now to transfer the approximation from B^ to B^, obtaining 
an approximation u^, while conserving the first integral, i.e. T[u^] = I[u^]. 
This can be formulated as a constrained minimization problem: 


mm 

uh^Bh 


\U 


-U^\\l2 S.t. Z\u^]=Z[u'^]. 


We observe that 


M M 


M M 


M M 


^ I Il2 — UiUjAij 2 UiUj Cij -|- ttj Uj Aij 

i=0 j=0 i=0 j=Q i=0 i=0 

= u^iu - 2u^Cu” + u^Au”, 
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where Aij = {ipi,ipj)i 2 , Aij = {(pi,(pj)L 2 and Qj = {(pi,(pj)L 2 . Also observ¬ 
ing that 


I[u'^] = Ipn + l{u), I[u^] = Xpn(u”), 

the problem can be reformulated as 

min u'^Au — 2u^Cu’^ -I- u”Au"' s.t. Xpn+i(u) — Xpn(u"') = 0. 
ueR"+i 

This is a quadratic minimization problem with one nonlinear equality con¬ 
straint. Using the method of Lagrange multipliers, we find u as the solution 
of the nonlinear system of equations 

iu - Cu^ - AVXpn+i (u) = 0 
Xpn + l(u) -Xpn(u"') = 0, 

which can be solved numerically using a suitable nonlinear solver. 

In general, applicable also in the FDM case, given u obtained by inter¬ 
polating u” onto in a non-preserving manner, a preserving transfer 

operation is obtained by solving the system of equations 

u — u — AVXpn+i (u) = 0 

Xpn+i(u) -Xpn(u") = 0. 


3.2 Projection methods 

Let the function /p : x be such that 


un+l _ 


= /pK,u-+i) 


At 


(3.2) 


defines a step from time tn to time tn+i of any one-step method applied 


to (2.1) on the fixed grid represented by the discretization parameters p. 


Then we define one step of an integral preserving linear projection method 
to p"+^ by 


u"" I—>■ from p" 


1. Interpolate u"' onto by whatever means to get u. 


2. Integrate u one time step by computing u = u -|- At/pn+i (u, u), 

3. Compute by solving the system of M -|- 1 equations = 

u -I- Az and Xpn+i = Xpn(u"'), for G and A G M, 

where the direction of projection z is typically an approximation to 
VXpn + l(u’"+l). 
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By utilizing the fact that for a method dehned by (|3.2|) there exists an 

jM 

P . ii’^ 7 iL\i. LilCXt IX - 


implicitly dehned map ^'p : —)■ such that = 'I'pU”, we dehne 


ffp(u"):= 


^'pU"^ — 

At ’ 


and may then write the tree points above in an equivalent, more compact 
form as: Compute G and A G M such that 


u 


n+l 


— u — Af^pu+i (u) — Az = 0, 

p 


s/p" 

Xpn + l(u"+l) -Xpn(u") = 0, 


(3.3) 

(3.4) 


where u is u"" interpolated onto by an arbitrary procedure. 

The following theorem and proof are reminiscent of Theorem 2 and its 
proof in 22 , whose subsequent corollary shows how linear projection meth¬ 
ods for solving ODEs are a subset of discrete gradient methods. 

Theorem 1. Let gp : —)■ be a consistent discrete approximation of f 

in (2.1) and /et VXp(u”, u”"*"^) be any discrete gradient of the consistent ap¬ 
proximation Ip{u) ofI[u] defined by (2.2) on the grid given by discretization 
parameters p. If we set S'pn+i in (3.1) to be 


Spn+l (u, U 


„+i, _ gp..+i(u)z^ - Zgp»+l(u)- 


(VXp 


1+1 (U, u 


n+l 


),z) 


(3.5) 


then the linear projection method for solving PDFs on a moving grid, given 
by is equivalent to the discrete gradient method on moving grids, 

as given by (3.1). 

Proof. For better readability, we set VX := VXpn+i (u, u”"*"^). Assume that 
(3.3)-(3.4) are satished. By applying (3.4), we get that 

Xpn(u’^) -Xp„ + l(u) =Xp„ + l(u”+l) -Xpn + l(u) 

= (VX, u”+^ - u) 

= At (VX, ^pn+i (u)) A (VX, z) , 


and hence 


^ ^ Xpn(u’") -Xpn + l(u) _ ^^ (VX, gp n + l(u)) 


(VX,z) 

Substituting this into ( |3.3[ ), we get 

Xpn(u’^) -Xp„+i(u) 


VX,z 


(3.6) 


= u + 


(VX,z) 


z -I- At 5 pu+i (u) - 


(VX,ffpn+l(u)) 

(VX,z) 
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where 


(VX,ffpn + l(u)) 


z = 


VX'^zg'pu+i (u) — VX'^( 5 rpn+i (u)z 
(VX,z) 

( 7 pn+i (u)z'^VX — zg^n+i (u)'^VX 
(VX,z) 


and thus (3.1) is satisfied, with S'pu+i as given by (3.5). Conversely, if 


satisfies (3.1), then (3.4) is satisfied. Furthermore, inserting (3.5) into (3.1) 


and following the above deduction backwards, we get (3.3), with A defined 
by (pi. □ 


the linear pro¬ 
jection methods on moving grids constitute a subset of all possible discrete 


Since (3.5) defines a particular set of choices for 5pn+i, 


gradient methods on moving grids as defined by (3.1). Note also that, since 


the linear projection methods are independent of the discrete gradient, each 


linear projection method defines an equivalence class of the methods (3.1), 
uniquely defined by the choice of ^pn+i. 


3.3 Family of discretized integrals 

At the core of the methods considered here is the notion that an approxi¬ 
mation to the first integral X is preserved, and that this approximation is 
dependent on the discretization parameters which may change from iteration 
to iteration. That is, we have a family of discretized first integrals Xp, and 
at each time step the discretized first integral is exchanged for another. For 
each set of discretization parameters p, there is a corresponding set of de¬ 
grees of freedom u, in which we search for a u such that Xp(u) is preserved. 
This can be interpreted as a fiber bundle with base space B as the set of all 
possible discretization parameters p, and fibers Xp as the sets of all degrees 
of freedom such that the discretized first integral is equal to the initial dis¬ 
cretized first integral, i.e. Xp = {u G M^|Xp(u) = Xpo(u*^)}. A similar idea, 
although without energy preservation, has been discussed by Bauer, Joshi 
and Modin in (2^ . 


4 Numerical experiments 

4.1 General remarks on type of experiments made 

To provide examples of the application of our method and to investigate 
its accuracy, we have applied it to two one-dimensional PDFs: the sine- 
Gordon equation and the Korteweg-de Vries (KdV) equation. The choice 
of these equations were made because they both possess traveling wave so¬ 
lutions in the form of solitons, providing an ideal situation for r-adaptivity. 
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which allows the grid points to cluster around wave fronts. The following 
experiments consider r-adaptivity only, and not p- or /i-adaptivity. The sine- 
Gordon equation is solved using the FDM formulation of section 2.2 while 
the KdV equation is solved using the PUM formulation of section [2^ 

We wish to compare our methods to standard methods on fixed and adap¬ 
tive meshes. This gives us four methods to consider: Fixed mesh methods 
with energy preservation by discrete gradients (DG), adaptive mesh meth¬ 
ods with preservation by discrete gradients (DGMM), a non-preserving fixed 
grid method (MP), and the same method with adaptive mesh (MPMM). The 
former two methods are those described earlier in the paper, while the latter 
two are made differently for the two equations. In the sine-Gordon case, 
we use a finite difference scheme where spatial discretization is done using 
central finite differences and time discretization using the implicit midpoint 
rule. In the KdV case, the spatial discretization is performed the same way 
as for the discrete gradient schemes, while the time discretization is done 
using the implicit midpoint rule. Note that the MPMM scheme for the sine- 
Gordon equation proved unstable unless restrictively short time steps were 
used, and the results of those tests are therefore omitted from the following 
discussion. The results using MPMM for the KdV equation were better, and 
are presented. The procedure for mesh adaptivity is presented in the next 
subsection. 

Using adaptive mesh methods, one obtains a different set of mesh points 
at each time step, meaning the numerical solution u from the previous time 
step must be transferred onto the new set of mesh points. We tested three 
different ways of doing this, two of which are using linear interpolation 
and cubic interpolation. The linear interpolation consists of constructing 
a function u(x) which is piecewise linear on each interval [xf, such that 
u{x^) = u2, then evaluating this function at the new mesh points, giving the 
interpolated values Ui = u{x^^^). The cubic interpolation consists of a simi¬ 
lar construction, using cubic Hermite splines through the MATLAB function 
pchip. Of these two transfer methods, the cubic interpolation yielded su¬ 
perior results in all cases, and so only results using cubic interpolation are 
presented. The third way, using preserving transfer as presented in section 
3.1, applies to the KdV example, where the PUM is used. Here, we found 


little difference between cubic interpolation and exact transfer, so results are 
presented using cubic interpolation for the transfer operation here as well. 


4.2 Adaptivity 

Goncerning adaptivity of the mesh, we used a simple method for r-adaptivity 
which can be applied to both FDM and FEM problems in one spatial di¬ 
mension. When applying moving mesh methods, one can either couple the 
evolution of the mesh with the PDE to be solved through a Moving Mesh 
PDE |24| or use the rezoning approach, where function values and grid points 
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are calculated in an intermittent fashion. Since our method is based on hav¬ 
ing a new set of grid points at each time step, and not coupling the evolution 
of the mesh to the PDE, the latter approach was used. It is based on an 
equidistribution principle, meaning that when 0 = [a, b] is split into M 
intervals, one requires that 


^i+1 


uj{x)dx ~ j ^{x)dx, 


Xi 


where the monitor function w is a function measuring how densely grid points 
should he, based on the value of u. The choice of monitor function is problem 
dependent, and choosing it optimally may require considerable research. A 
variety of monitor functions have been studied for certain classes of problems. 


see e.g. 25 26 . Through numerical experiments, we found little difference in 


performance when choosing between monitor functions based on arc-length 
and curvature, and have in the following used the former, that is, the gener¬ 
alized arc-length monitor function |25] 


In this case, the equidistribution principle amounts to requiring that the 
weighted arc length (in the case k = 1 one recovers the usual arc length) of u 
over each interval is equal. In applications, we only have an approximation 
of u, meaning uj must be approximated as well; in our case, we have applied a 
hnite difference approximation and obtained approximately equidistributing 
grids using de Boor’s method as explained in |15[ pp. 36-38]. We tried 
different smoothing techniques, including a direct smoothing of the monitor 
function and an iterative procedure for the regridding by De Boor’s method 
(see e.g. 27 ^). In the case of the KdV equation, there was little to 


no improvement using smoothing, but the sine-Gordon experiments showed 
signihcant improvement with direct smoothing; i.e., in De Boor’s algorithm, 
we use the smoothed discretized monitor function 

OJi-l + 2LOi + 

UJi = -^-. 


4.3 Sine-Gordon equation 

The sine-Gordon equation is a nonlinear hyperbolic PDE in one spatial and 
one temporal dimension exhibiting soliton solutions, with applications in 
predicting dislocations in crystals and propagation of fluxons in junctions 
between superconductors. It is stated in initial value problem form as: 

utt-Uxx + sm{u) = {x,t) gRx [0,T], (4.1) 

uix,0) = f{x), utix,0) = g{x). 
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We consider a finite domain [—L, L] x [0, T] with periodic bonndary condi¬ 
tions u{—L) = u{L) and Ut{—L) = Ut{L). The eqnation has the hrst integral 


I[u] 


1 2 1 2 . 

— 'lif “h 1 

2 c ' 2 ^ 


cos(n)dx. 


Introdncing v 


Ut, (4.1) can be rewritten as a hrst-order system of PDEs: 


Ut 


V 



Uxx - sin(?x) 


with hrst integral 


I[u, v] 


-v^ H —ui -|- 1 — cos(u)dx. 
2 2 


(4.2) 


Finding t he v ariational derivative of this, one can interpret the eqnation in 
the form (2.3) with S and as follows: 


S = 


We will apply the FDM approach presented in section |2.2[ approximating 
(4.2) by some qnadratnre with points and weights 


■ 0 

1' 

5X 

, — \u, v\ = 

sin(u) - Uxx 

-1 

0 


V 


M 


I[u, n] ~ ^ i 


i=0 


+ 1 - cos{ui) 


In addition, we approximate the spatial derivatives with central differences. 
At the endpoints, a periodic extension is assnmed, yielding the approxima¬ 
tion 


Here, 6wi = tCj+i — Wi-i denotes central difference, with special cases 6uo = 
6 um = ui — um- 1 , and 5xo = Sxm = xi — xq + xm — xm-i- Taking the 
gradient of Tp(u) and applying the AVF discrete gradient gives 


1 

VXp(u^, u”+i) = j VXp(^u” + (1 - 
0 

The periodic bonndary conditions are enforced by setting uq = um- In the 
implementation, the k* were chosen as the qnadratnre weights associated 
with the composite trapezoidal rnle, i.e. 


Kq = 


xi - Xq 


Km = 


Xm — xm-1 


Ki = 


Xi-\-l Xi—i 


i = 1, ...,M - 1. 
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Figure 1: Left: Illustration of kink-antikink solution. Right: Grid movement 
- each line represents the path of one grid point in time. 


Furthermore, S was approximated by the matrix 


Sd 


0 I 
-I 0 ’ 


with I an M X M identity matrix. The exact solution considered was 


u{x,t) = 4 tan ^ 


/ 

V 


sinh f ] 

c cosh ( , I 


\ 


This is a kink-antikink system, an interaction between two solitons, each 
moving in different directions with speed c G (0,1), resulting in two wave 
fronts traveling in opposite directions. The wave fronts become steeper as 
c —)> 1. 

Figure [^illustrates the analytical solution and shows the time evolution 
of the mesh as obtained with the DGMM method. Note that the grid points 
cluster along the wave fronts. The left hand side of Figure shows the 
time evolution of the error = ||u^(x) — u{x,tn)\\L 2 ^ where is a linear 
interpolant created from the pairs (u”,x"'). The right hand side of Figure 
[^ shows the time evolution of the relative error in the discretized energy, 
F;^ = (V(u-)-/po(uO))//po(uO ). We can see that the long-term behaviour 
of the MP scheme is superior to that of the DG scheme, but when mesh 
adaptivity is applied, the DGMM scheme is clearly better. Also note that 
while the DG and DGMM schemes preserve Ip to machine precision, the MP 
scheme does not. 

Figure [^ shows the convergence behaviour of the three schemes with 
respect to the number of spatial discretization points M, and the number of 
time steps N. Note that the DG and MP methods plateau at ~ 400; this 


18 






















































































t 



Figure 2: Left: L 2 error. Right: Relative error in/p. Parameters: At = 0.01, 
M = 300, L = 30, c = 0.99. 



Figure 3: Left: Error at T = 8 as a function of M, with At = 0.008, c = 0.99, 
L = 30. Right: Error at T = 8 as a function of iV = T/At, with M = 1000, 
c = 0.99, L = 30. 




epsilon 

Figure 4: Error at T = 8 as a function of e, with At = 0.01, M = 600 and 
L = 30. 
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is likely due to the error stemming from spatial discretization dominating 
the time discretization error for these methods, while the DGMM scheme 
has lower spatial discretization error. 

Finally, to illustrate the applicability of the DGMM scheme to harder 
problems. Figure shows the error at stopping time of the methods as a 
function of a parameter e representing the increasing speed of the solitons 
(c = 1 — e). From this plot, it is appararent that while the non-adaptive 
MP scheme is competitive at low speeds, the moving mesh method provides 
significantly more accuracy as c —>■ 1. 


4.4 Korteweg—de Vries equation 

The KdV equation is a nonlinear PDF with soliton solutions modelling shal¬ 
low water surfaces, stated as 

Ut Uxxx T QuUx — 0* (^■^) 


It has infinitely many first integrals, one of which is the Hamiltonian 



— u^dx. 


With this Hamiltonian, we can write (4.3) in the form (2.3) with S and 
as follows: 


m 

Su 


s = 


d_ 5n 

dx’ 6u 


[w] = Uxx + 3u^. 


We will apply the PUM approach to create a numerical scheme which pre¬ 
serves an approximation to 'H[u], splitting H = [—L,L] into M elements 
{[xj, Xj_|_i]}^g ^ and using Lagrangian basis functions (pj of arbitrary degree 


for the trial space. Approximating u by as in section 

Wp(u) = Min'-] = 


2.3 


we find 


^ J,k,l 


(pj(pk(pidx. (4.4) 


The integrals can be evaluated exactly and efficiently by considering elemen¬ 
twise which basis functions are supported on the element before applying 
Gaussian quadrature to obtain exact evaluations of the polynomial integrals. 
We define 


Dijk — 


ipjipkifidx and Eij = / ipj^x<Pk,xdx. 
; Jq 
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The matrices A and B with 


Aij — 


(picpjdx and 


Bji — 


In 




are formed in the same manner. Note that B is in this case independent of 
u. Applying the AVF method yields the discrete gradient 


ViWp(u^, u”+i) = J + (1 - 
0 


such that, with the convention of summation over repeated indices, 


(V«p), = Iem + - d.mK + ‘ 


(§<+-r'))- 


This gives us all the required terms for forming the system (2.15) and apply¬ 
ing the discrete gradient method to it. During testing, the ipj were chosen 
as piecewise linear polynomials. The exact solution considered is of the form 


u{x, t) 


|sech2 , 


(4.5) 


which is a right-moving soliton with c as the propagation speed, chosen 
as c = 6 in the numerical tests. We have considered periodic boundary 
conditions on a domain [—L,L] x [0, T], with L = 100 in all the following 
results. 

Our discrete gradient method on a moving mesh (DGMM) is compared 
to the same method on a static, equidistributed mesh (DG), and the implicit 
midpoint method on static (MP) and moving mesh (MPMM). The spatial 
discretization is performed the same way in all cases. Figure shows an 
example of exact and numerical solutions at t = 15. Note that the peak in 
the exact solution will be located at x = ct. 

To evaluate the numerical solution, it is reasonable to look at the distance 
error 


= ctn - X*, 


where x* = argmaxu/i(x, he. the location of the peak in the numerical 

X 

solution. Another measure of the error is the shape error 



Uh{x,tn) - u 



where the peak of the exact solution is translated to match the peak of the 
numerical solution, and the shapes of the solitons are compared. 
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Figure 5: Solutions at T = 15. At = 0.01, M = 400. MP and DG are almost 
indistinguishable. 
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Figure 6: Relative error in the Hamiltonian plotted as a function of time 
t E [0,15]. At = 0.01, M = 400. 


Figure confirms that the DG and DGMM methods preserve the ap¬ 
proximated Hamiltonian (4.4), while it is also worth noting that in the case 
of the midpoint method, the error in this conserved quantity is much larger 
on a moving than on a static mesh. Similar behaviour is also observed for 
a moving-mesh method for the regularized long wave equation in the recent 
paper 29 , where it is concluded that a moving mesh method with a conser¬ 
vative property would be an interesting research topic. Figure where the 
phase and shape errors are plotted up to T = 15, is an example of how the 
DGMM method performs comparatively better with increasing time. 


In figures and we present the phase and shape errors for the differ¬ 
ent methods as a function of the number of elements M and the number of 
time steps N, respectively. We observe that the DGMM scheme performs 
especially well, compared to the other three schemes, for a coarse spatial 
discretization compared to the discretization in time. In figure the phase 
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Figure 7: Phase error (left) and shape error (right) as a function of time. 
At = 0.01, M = 400. 
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Figure 8: Phase error (left) and shape error (right) as a function of the 
number of elements M, at time T = 5. At = 0.01. 


and shape errors are plotted as a function of the parameter c in the exact 
solution (4.5), where we note that | is the height of the wave; increasing c 
leads to sharper peaks and thus a harder numerical problem. As expected, 
the advantages of the DGMM method is less evident for small c, but we ob¬ 
serve that the DGMM method outperforms the static grid midpoint method 
already when c = 2. 


4.5 Execution time 

The code used is not optimized, so any quantitative comparison to standard 
methods has not been performed; it is still possible to make some qualitative 
observations. Adding adaptivity increases time per iteration slightly since 
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Figure 9: Phase error (left) and shape error (right) at time T = 5, as a 
function of the number of time steps N = T/At. M = 800. 



Figure 10: Phase error (left) and shape error (right) as a function of c in the 
exact solution (4.5), at time t = 5. At = 0.01, M = 800. 
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the systems become more complicated, especially in the case of the PUM ap¬ 
proach where the matrices A and B need to be recalculated, at each time step 
when adaptivity is used. This increases runtime somewhat when compared 
to hxed grid methods. However, adaptivity allows for using fewer degrees 
of freedom, and so decreases the degrees of freedom needed for a given level 
of accuracy. This accuracy gain is more pronounced the harder the problem 
is (steeper wave fronts etc.), and so it stands to reason that there will be 
situations where adaptive energy preserving methods will outperform non- 
adaptive and/or non-preserving methods. This is in accordance with what 
we have observed from our not optimized experiments. 


5 Conclusion 


In this paper, we have introduced a general framework for producing adap¬ 
tive first integral preserving methods for partial differential equations. This 
is done by hrst providing two means of producing hrst integral preserving 
methods on arbitrary hxed grids, then showing how to extend these methods 
to allow for adaptivity while preserving the hrst integral. Numerical testing 
shows that moving mesh methods coupled with discrete gradient methods 
provide good solvers for the sine-Gordon and Korteweg-de Vries equations. 
It would be of interest to apply the method to higher-dimensional PDEs 
with a more challenging geometry, preferably using the PUM approach, to 
investigate its accuracy as compared to conventional methods, and to test 
whether h- and/or p-rehnement provides a notable improvement. It may 
also prove fruitful to explore the ideas presented in |23] to make the trans¬ 
fer operations between sets of discretization parameters in a more natural 
setting than simply interpolating, as suggested in section 3.3 Furthermore, 


analysis of the methods considered here could provide important insight into 
e.g. stability, consistency and convergence order. 
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